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Kinetic models described by systems of linear differential equations can be fitted to data quickly and easily by 
taking advantage of the special properties of such systems. The estimation situation can be greatly improved 
when multiresponse data are available, since one can then automatically determine starting values and better 
discriminate between rival models. 
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1. Introduction 

In this article we summarize the work of a series of 
papers [1-3] 1 in which we deal with Fitting First order 
kinetic models to uniresponse and multiresponse data. 
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of Mathematics and Statistics, Queen's University at 
Kingston. The work reported on was supported by the 
National Science Foundation in the United States and 
the Natural Sciences and Engineering Research Council 
of Canada. 

1 Numbers in brackets indicate literature references. 



We consider systems in which the expected responses 
at K points in the system, ij(O=0>ji(O, ""lit'). ■■• ^(0)'. 
are described by the system of linear differential equa- 
tions 



di}/dt=i)(t)=A i)(t)+i(t) 



0-1) 



where A is a K XK system transfer matrix depending on 
rate constants t , d it ... , and i(f) is a vector input function 
to the system. We assume further that there are K initial 
conditions ijo = (tjoi, no 2 , ■•■ , iio K y, some possibly un- 
known, and that r=t—9o, where 0o is a (possibly un- 
known) time delay. All the unknown parameters are 
gathered into aPxl parameter vector 0. 
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Example: Oil Shale. 1 

As aa exacple of a ehemica! system described by 3 set 
of linear differentia] equations, we cite the pyroiysis of 
oil shale in which the model, Fitted by Zlegel and 
Gorman [4] has the system diagram 



6 4 



\ 



n 3 



In this system, 77, denotes kerogen, % bitumen, and 173 
oil. The model implies that kerogen decomposes to bitu- 
men, with rate constant Jf aad. to oil vrilh rate constant 
9 4 , and bitumen produces oil with rate constant 2 > and 
unmeasured by-products with rate constant 6 3 . 



2. Expectation Functions and Derivatives 

With Respect to the Parameters for First Order 

Kinetic Systems 

Several methods are used to estimate the parameters 
in first order kinetic models. The most obvious method 
is to solve the system of differential equations corre- 
sponding to the particular compartment model and use 
the resulting expectation function in a standard non- 
linear estimation program, A second approach is to fit a 
genera] suei of exponentials model by "peeling" [51 A 
Ihirn approach is to use a standard nonlinear estjmaii&n 
program, using numerical integration to solve the equa- 
tions. A superior approach, proposed by Jennrich and 
Bright [61, is to obtain the general solution to the system 
of equations by calculating values for the model func- 
tion tj(0 and its derivatives directly, given values of 9 
and i and i(f). 

2,1 The General Solution 

The solution to a linear system of diffeienti&i equa- 
tions can be expressed in terms of convolutions using the 
matrix exponential (7J. The solution is 



where the **' denotes convolution, 



e- J "i(0=V-^(3d£. 



(2-1) 



{1.1) 



The integral of the vector function is evaluated com- 
ponentwise, Gampsitariooal methods fat ■equating the 
convolution integral are gh'&n in Moier and Van Loan 
[8], when A is diagonahzable, and in Bavely aad Stewart 

[9], when A is nondiagona3izaWe. 



2.2 Derivatives <rf the Expectation Function 

To use a Gauss-Newton procedure to estimate the 
parameters we need derivatives with respect to the pa- 
rameters. As shown by Jsnnrich and Bright [6], a. great 
advantage to compartment models is that the deriva- 
tives can be evaluated jn the same fashion as the model 
function itself. Instead of differentiating 3j(r) directly as 
in Jennrich and Bright, however, we differentiate eq 
(1.1) and solve the resulting linear system of differential 
equations. This idea was discussed in another context in 
Smith [10! and was ^^ by KEilbfLeisch et ah [11]. 

To simplify notation, we use a subscript p to denote 
differentiation with respect to the parameter b ? 30 



%to= 



318(0 



»0, 



for p = 1, X ... P. The derivative of (1,1) with respect to 



Q p is then 



fc(0 ^Att,{f) +A f i){t)+t,{t} 



for which the solution is 

y p (t) = e A \(V)+e M *lA,r}(t) + t f (m (2.3) 

Dead time, a , can he incorporated by modifying ( 1 . 1) to 
r[(r)=Av\ij}+t 

■V(0) = TJO 



where 



r=s 



r— 0o i>Qo 







t<0o. 



l£@o is known, we simply replace f by r in eqs 1.1,2.1, 
and 2.3, but if So is unknown and depends on a parame- 
ter, the expression for the derivatives is extended to 



V P (r)^e^ Vp (0ne Ar *[A^(7)+ h ] 
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It is easy to evaluate A p , i\ p (0) = a ij o/a $ p , and t p since 
they are constants and, in fact, usually two of the three 
are zero. Also, for any t<do, t p is zero for all 
p = l, ... » P. Note that the method can be extended to 
higher order derivatives. 

Specifying the Model 

A simple unambiguous computer notation can be used 
to specify first order kinetic models in & parameter table 
consisting of three columns, the first column giving the 
parameter number. For a rate constant, the second col- 
umn entry gives its source and the third column entry its 
sink, a sink with compartment number denoting elimi- 
nation. For initial conditions, tfo, the second column 
entry is the number of the component in i}o and the third 
column entry is — 1. For a step input, i, the second 
column entry is the number of the component in i and 
the third column entry is —2. Dead time is coded as in 
column 2 and in column 3. 

Example: Oil Shale.2 

The oil shale parameter table is presented in annotated 
form. It is to be noted that a single parameter may repre- 
sent more than one rate constant. 

Table 1. Oil shale parameters. 



Parameter 
number (type) 



Column 



Column 



1 


(rate constant) 


1 


(source) 


2 


(sink) 


2 


(rate constant) 


2 


(source) 


3 


(sink) 


3 


(rate constant) 


2 


(source) 





(elimination) 


4 


(rate constant) 


1 


(source) 


3 


(sink) 


5 


(dead time) 













3. Multiresponse Estimation 

In the multiresponse situation when the errors have 
unknown variances and covariances but are assumed to 
be temporally uncorrelated, the appropriate criterion 
derived via a likelihood or Bayesian approach is to min- 
imize the MxM determinant [12]. 



V(0)\ = \Z'Z\ 



(3.1) 



In eq.(3.1), Z = Y-H is the NxM matrix of residuals 
{z n J={v*0«)-Ut ('-.Xh k = \, ,2 , ... , M, n=\, 2, ... , 
N. The expected responses ij are assumed to depend on P 
parameters 6; except where explicitly required, we sup- 
press the dependence on 6. For first order kinetic sys- 
tems with all responses measured, M=K, 



To determine the minimum of | Z'Z \ , it is advan- 
tageous to calculate the gradient and Hessian and ex- 
ploit Gauss-Newton optimization techniques. Efficient 
numerical procedures for computing the gradient and 
approximate Hessian are given in [2], and an algorithm 
which performs the calculations in [1], 

Alternative expressions for the components of the 
gradient y and the Hessian T are 



y p = d\V\/dd p = 2\y\tr[V- [ Z'Z p ] , 
p = l,2,...,P , 



and 



r w = e 2 |F|/afV0,= 



(3.2) 



(3.3) 



= y p y,+2 | VltriV'Z'Z^-'Z'Z,] 
+2 1 V \tr[V- x Z'Z,V~ l Z;Z] 

+2\V\tr[V- 1 Z p Z q ] + 2\V\tr[V- l Z'Z pii ) , 
p,q = \,2,...,P . 

The second derivative terms Z n in eq (3.3) are ignored 
to produce an approximate Hessian. 



4. Practical Aspects 

Linear Constraints 

Sometimes the data matrix Y involves dependencies 
as a result of imputation of responses or mass-balance 
calculations. If these dependencies also occur in the ex- 
pected responses, then important modifications to the 
multiresponse estimation procedure must be made so as 
to avoid convergence to spurious optima [13,14]. It is 
therefore necessary to examine the residual matrix Z(6) 
for singularities, which can be done by arranging the 
rounding units in the columns of Y to be approximately 
equal and taking a singular value decomposition of Z 
[15]. As explained there, singular values on the order of 
the rounding unit indicate singularity and should 
prompt the analyst to search for constraints in the data. 
Such examinations should be done at the beginning of 
the analysis using the initial parameter values and at the 
end of the analysis using the converged values. To aid 
convergence, logarithms of parameters are used during 
estimation. 

Linear constraints can be dealt with easily by com- 
bining the linear constraint vectors into a matrix, per- 
forming a QR decomposition of that matrix, and letting 
the rotation matrix W be the columns of Q which are 
orthogonal to the constraint vectors. We then simply 
minimize \{ZW)\ZW)\, where ZW=YW-HW. 
Clearly, the gradient and Hessian of this determinant 
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are obtained from eqs 3.2 and 3.3 by replacing Z by ZW 
and Z p by Z„W. 

Constraints on the Number of Parameters, 
Responses and Observations 

The determinant criterion implies two constraints on 
the number of observations [2,3]. First, N must be at 
least equal to M since otherwise the determinant is iden- 
tically zero. Second, N must exceed P otherwise the 
criterion can be made zero by fitting any one response 
perfectly, which can generate up to M distinct minima. 
Thus the residual matrix has effectively TV — P degree of 
freedom. It may seem that there should be more degrees 
of freedom since there are NM separate observations, 
but the criterion can be locally controlled by any one 
response so the effective number of observations is N 
rather than NM. 

Starting Values 

An important part of fitting nonlinear models is deter- 
mining good starting values. For uniresponse data, an 
effective method is to use peeling in which we plot the 
logarithm of the response versus time and fit a straight 
line to the segment at large t values. The slope of the line 
gives an estimate of the smallest eigenvalue of the A 
matrix. Using the fitted line to generate residuals and 
plotting the logarithm of the residuals versus t should 
again reveal a straight line portion at large t values, so 
the process is repeated, thereby obtaining estimates for 
the eigenvalues. As mentioned in section 2, this process 
is often used for parameter estimation, but we do not 
recommend it. 

In the case of multiresponse data for first order kinet- 
ics, the problem is easily solved using linear least squares 
by exploiting the linear relation between the rates and 
the responses! As noted in [3], if we could measure the 
rates y and the responses y at a particular time r, then 
using y(r) =A y(r) produces a linear relation between 
the "dependent" variable y=-y and the "independent" 
variables x p =A p y in the form y =X6. We can thus 
solve for 6 by using linear least squares. A simple pro- 
cedure for obtaining starting values, then, is to use ap- 
proximate rates from Finite differences of the responses 
at successive time points and x p values from the corre- 
sponding averages. Alternatively, one could smooth the 
data for each response by fitting splines so as to obtain 
better rate and response values, and then use these in a 
linear least squares routine. 

Example: a-pinene.l 

Data on the thermal isomerization of a-pinene at 
189.5° were reported by Fuguitt and Hawkins [16], and 



a thorough multi-response analysis was presented in 
[13]. The fitted model was described by the system 



7 = 



which can also be written y=X8, where 



7i 




-0i7i 


-0 2 y. 






72 




0,7. 








73 


= 




2 7i- 


-0373-0473+0575 


= Ay 


74 








0373 




75 








0+73 -0s7s 





X = 



r~-7i -7i 

7, 

0, 73 











-7i -7s 

73 

7 - 






y 
o 

■7s 



Substituting estimated rates for each time and joining 
them into a single vector y, and calculating X matrices 
for each time and joining them into a single X matrix, 
allows us to use linear regression to estimate starting 
values for 0. The starting values obtained are listed in 
column 2, table 2. Note their closeness to the converged 
values, column 4. 

Table 2. Estimation for a-pinene at 189.5°. 







Box 


Bates and Watts 


Parameter 


Start 


etal. 


Full 
Model 


Reduced 
Model 


1 


5.84 


5.95 


5.94 


5.95 


2 


2.65 


2.85 


2.84 


2.82 


3 


1.63 


0.50 


0.45 


- 


4 


27.77 


31.5 


31.21 


30.75 


5 


4.61 


5.89 


5.79 


5.72 


Determinant 












600 




28.4 


29.0 



Parameter estimates are listed in table 2 in 
minutes" 1 X 10 _s ; that is a table value of 5.84 is actually 
a rate constant of 5.84X 10 -5 minutes -1 . 

When some of the responses are not measured, it is 
still possible to use approximate rates provided other 
information, such as a mass-balance, is substituted. 



5. Two Examples 



5.1 Oil Shale 



The model and data presented by Ziegel and Gorman 
[4] were fitted using the procedures described here. In 
this problem the concentration of tj, was not measured, 
which introduces some complexity in determining start- 
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ing values. Second, the elimination compartment, corre- 
sponding to coal and gas, was not measured. Third, 
there was a time delay caused by the shale having to 
reach the reaction temperature. The observed responses 
are>»2 and yi , measured in percent of the initial kerogen 
7), , so K = 3, M = 2, with N = 14. 

To determine starting values in this case, we plotted 
the data and obtained a rough estimate of 0o=6 min. 
Estimating the initial slopes of tjj and t) 3 from the graph 
gave values of 0.029 min -1 for 6 Y and 0.013 min -1 for + . 
This delay estimate and the rate constants were used to 
obtain starting estimates of 6 2 and 3 as shown in column 
2, table 3, following the procedure of section 4. The final 
results, together with those of Ziegel and Gorman, are 
shown in columns 3 and 4 of table 3. 

Table 3. Parameter estimates for oil shale data. 



Parameter 


Start 


Bates-Watts 


Ziegel-Gorman 


1 


0.029 


0.0172 


0.0173 (ki) 


2 


0.012 


0.0090 


0.0092 (kjfi) 


3 


0.028 


0.0205 


0.0201 <k 2 (t-f 2 )) 


4 


0.013 


0.0104 


0.0104 (k 3 ) 


00 


6.0 


7.8 


7.7 


Determinant 










67980 


429 





5.2 a-pinene 

In their analysis of the a-pinene data, Box et al. [13] 
noted that response 4 was imputed using 
,y 4 =0.03(1 00 -jr t ), and that the data set was subject to a 
mass-balance constraint, y]+y 2 +yi+y4+y5=100. To 
avoid convergence to spurious optimal parameter val- 
ues, they recommended that these data dependencies be 
taken into account by using observation vectors consis- 
ting of linear combinations of y t ,y lt y^,y\ andj>5 which 
are orthogonal to the space defined by the vectors 
(0.03, 0, 0, 1, 0)' and (1, 1, 1, 1, 1)'. We therefore 
treaty as an unmeasured component for estimation pur- 
poses, and use linear combinations of the responses 
which are orthogonal to the vectors al = (0, 0, 0, 1, 0) 
and a 2 =(l, 1, 1, 1, 1). The rotation matrix M and the 
modified responses can therefore be determined by per- 
forming a QR decomposition on the matrix (a,, a 2 ) and 
using the last 3 columns of Q coupled with all the re- 
sponses. In this case, K = 5, M = 3, with TV = 8. 

Approximate 95% confidence limits for 1«0 3 were 
very wide, suggesting that 9 3 was badly estimated and 
could be zero. We therefore fitted a reduced model 
in which there was no path from tj 3 to i) 4 , see column 5. 
The change in the determinant is 0.6 on 1 degree of 
freedom, which, when compared with the scaling factor 
s 2 =28. 4/3— 9.46 on 3 degrees of freedom, is clearly 



small, verifying that the reduced model is adequate for 
this data set. 

To further substantiate the adequacy of the reduced 
model, we fitted both models to a second set of data 
taken 204.5° [16]. The results of this fitting procedure 
are presented in table 4. The reduced model appears to 
be adequate for both data sets. 

Table 4. Estimation for a-pinene at 204.5°. 



Parameter 



Start 



Full 
Model 



Reduced 

Model 



1 
2 
3 
4 
5 
Determinant 



23.0 
13.2 
8,3 
76.9 
16.0 

116 



22.6 
12.6 
0.02 
72.8 
15.6 

0.55 



22.6 
12.6 

72.8 
15.6 

0.55 



6. Conclusions 

Several advantages of the direct multiresponse esti- 
mation approach for systems of differential equations 
are apparent. First, the model can be specified directly 
from the network diagram. Second, there is no need to 
obtain the analytic solution to the differential equations 
describing the reactions. Third, there is no need to code 
the model functions in a nonlinear estimation routine. 
Fourth, the bothersome and error-prone step of obtain- 
ing and coding derivatives of the expected responses 
with respect to the parameters is eliminated. Fifth, ex- 
cellent starting values can be determined automatically. 



References 

[1] Bates, D. M., and D. G. Watts, A multi-response Gauss-Newton 
algorithm, Commun. Statist. — Simul. Comput, B(13), 
705-715 (1984). 

[2] Bates, D. M., and D. G. Watts, A generalized Gauss-Newton 
procedure for multi-response parameter estimation, SIAM 
Journal on Scientific and Statistical Computing, to appear. 

[3] Bates, D. M., and D. G. Watts, Multiresponse estimation with 
special application to systems of linear differential equations, 
Technometrics 27, 190-201 (1985). 

[4] Ziegel, E. R., and J. W. Gorman, Kinetic modelling with multi- 
response data, Technometrics 22, 139-151 (1980). 

[5] Anderson, D. H-, Compartmenta! modeling and tracer kinetics, 
Springer-Verlag (1983). 

[6] Jennrich, R. 1., and P. R. Bright, Fitting systems of linear differ- 
ential equations using computer generated exact derivatives, 
Technometrics 18, 385-392 (1976). 

[7] Noble, B., Applied linear algebra, Prentice Hall (1969). 

[8] Moler, C, and C. Van Loan, Nineteen dubious ways to compute 
the exponential of a matrix, SIAM Review 20 (1978). 



437 



[9] Bavely, C A., and G. W, Scew^rC, An a'.gorilbm for computing 
reducing subspaces by block diagonafizalion, S1AM Journal 
of Numerical Analysis 16, 359-367 (1979). 

[10] Smith, W. R., Parameter estimation in nonlinear models of bio- 
logical systems, Fisheries and Marine Services Technical 
Report S89, Fisheries and Environment Canada (1979). 

[1 1] Kalbfleisch , J. D.; J. F . Lawless and W. M. Vollmer, Estimation 
in Markov models from aggregate data, Biometrics 39, 
907-919 {\%i), 

[12] Box, G. E. P., and N. R. Draper, The Eayesian estimation of 
common parameters from several responses, Biomelrika 52, 
355-3S5 (1965). 



[13] Box, O. E. P.; W. G. Hurler, J. F MacOregcr, and J. Erjavec, 

Some problems associated with the analysis of raultiresponse 

data, Technometrica 15, 33-51 (1973), 
[14] McLean, D. D-; D. S. Pritchard, D. W. Bacon, and J. Downie, 

Singularities in multiresponse modelling, Technometrics 2% 

291-29& (1979), 
[IS] Dongarra, J. J.; J. R. Bunch, C. B. Moler, and G. W. Stewart, 

Liiipack User's Guide, Chap. II, Philadelphia-. S.LA.M. 

(1979). 
[i 6] Fuguut, R, E., and J. E. Hawkins, Rate of thermal isa merization 

of a-pineae in the liquid phase, J. Amer. Chem. Soc. 69, 

319-312(1947). 



DISCUSSION 

of the Bates- Watts paper, 
Multiresponse Estimation With Special 
Applications to First Order Kinetics. 

Michael Frenklach 
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The authors presented an interesting approach to pa- 
rameter estimation for first order kinetic systems. The 
method is user oriented and particularly suited for com- 
puter implementation as a "canned" program. Indeed, 
present chemical kinetic codes input reaction mech- 
anism in a natural chemical language, that is, specifying 
reactions (usually in unformatted READ routines) as 
they are conventionally written on the paper. This infor- 
mation is automatically converted to a so-catted reac- 
tion matrix and, based on it, to differential equations 
describing the kinetics of reaction species. The reaction 
matrix, which contains all the stoichiometry of the sys- 
tem, can conveniently provide the required input infor- 



1 Michael Frenkiach's contribution to the subject stems from work 
performed in the Department of Chemical Engineering, Louisiana 

State Univmtiy. 



Another important feature, from the user's point of 
view, is that the presented method is applicable to mul- 
tiresponse data. It should be realized that modern prob- 
lems of interest to chemical kinetics get tougher, as for 
example, formation of pollutants in hydrocarbon com- 
bustion. The experimental answer to the growing com- 
plexity of the systems is the employment of multiple 
diagnostics For simultaneous monitoring of various pro- 
cess variables. However, interpretation of the experi- 
mental results cannot be fully realized without reliable 
and convenient multiresponse methods. 

The following are some of my thoughts on the needs 
in this area: 

1) Often, kineticists exhibit a philosophical re- 
sistance to a multiparameter approach to experi- 
malion for automatic coding of the method of Bates and 
Watts. 
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